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Abstract 

Evolution is a dynannic process. The two classical forces of evolution are mutation and selection. As¬ 
suming small mutation rates, evolution can be predicted based solely on the fitness differences between 
phenotypes. Predicting an evolutionary process under varying mutation rates as well as varying fitness 
is still an open question. Experimental procedures, however, do include these complexities along with 
fluctuating population sizes and stochastic events such as extinctions. We investigate the mutational 
path probabilities of systems having epistatic effects on both fitness and mutation rates using a theo¬ 
retical and computational framework. In contrast to previous models, we do not limit ourselves to the 
typical strong selection, weak mutation (SSWM)-veg\rr[e or to fixed population sizes. Rather we allow 
epistatic interactions to also affect mutation rates. This can lead to qualitatively non-trivial dynamics. 
Pathways, that are negligible in the SSWM-regime, can overcome fitness valleys and become accessible. 
This finding has the potential to extend the traditional predictions based on the SSWM foundation and 
bring us closer to what is observed in experimental systems. 
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INTRODUCTION 


How repeatable is evolution? As the metaphor by Stephen J Gould goes ‘if we run the tape 
of life back from the start how likely is it that we will get the same outcome that we see around 
us today?’ [1], The pioneering work of Lenski et al. tackled this question experimentally with 
microbes. It is now possible to literally play back evolution from a certain starting point and see 
where it leads [2-6]. 

Such empirical explorations made the until then theoretical concept of fitness landscapes tan¬ 
gible. The concept of a fitness landscape is a mapping between the genotype and the phenotype 
of an organism. Since selection acts on the phenotype or essentially on the fitness of the phe¬ 
notype, the genotype of each phenotype can be attributed a certain fitness. Connecting the 
genotypes which are one mutational step away from each other leads to the concept of fitness 
landscapes [7, 8]. Such empirical studies do make it clear that predictions will not be based on 
simple rules but complicated phenomena such as epistasis and epigenetics which play a major 
role in the process of evolution [6, 9, 10]. 

Epistasis is any deviation from the additive effects of alleles at different loci [11]. Epistasis 
gives rise to rugged fitness landscapes which have been found to be quite common in experimental 
observations in a variety of model systems [12, 13]. In particular, reciprocal sign epistasis is a 
necessary condition for having a rugged fitness landscape [14]. While in magnitude epistasis the 
fitness always increases (or decreases) with every additional mutation in a non-additive manner, 
in sign epistasis, however, valleys appear in the fitness landscape. A certain mutation might 
have a lower fitness than the previous state although it leads to higher fitness eventually. In 
such a case not all paths in the fitness landscape might be accessible by the population [15]. 
Comparing experimental systems to theoretical predictions made on the basis of the underlying 
fitness landscape helps elucidate the role of microscopic properties of the system in determining 
the macroscopic evolutionary trajectory. The details of the process such as the mutation rate, 
fitnesses of individual states and the global population size act as constraints on the accessibility 
of paths [13]. Using the assumption of strong selection and weak mutation rates (SSWM), the 
system advances on the fitness landscape in a stepwise fashion. This automatically limits the 
possible number of adaptive paths [10]. 

Evolutionary predictability and the speed of the dynamics is not only determined by the molec¬ 
ular constraints of fitness and mutation rate but also by population dynamics [14]. Theoretical 
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explorations often assume a fixed population size starting at one node of the fitness landscape 
and its movement is tracked over the course of time. Increasing the population size, or the 
mutation rate, we observe the phenomenon of clonal interference [15, 16]. This occurs when 
a second step mutant arises in a population even when the first step mutation is not fixed. In 
other words, the SSWM assumption is no longer valid. Clonal interference has been extensively 
explored experimentally [17-19] as well as theoretically [16, 20-25]. This phenomenon removes 
the limit on the accessibility of non-adaptive trajectories. If the fitnesses and mutation rates 
align to particular conditions, i.e. the mutation rates also underlie epistatic interactions, then 
such valley crossings might be faster than adaptive trajectories [24, 26]. 

Populations in real systems are finite and their size can undergo fluctuations which can lead to 
possible extinction events. Together with the phenomena of clonal interference and epistatic in¬ 
teractions between mutations (correlated rugged fitness landscapes), predicting evolution through 
a given fitness landscape seems like an impossible task. Herein we develop a general methodol¬ 
ogy for predicting the most probable path in a fitness landscape with epistatic interactions in a 
multi-dimensional fitness landscape. To reflect a realistic scenario we use a multi-type branching 
process (e.g. [27]) to drop the assumption of a constant population size. For presentation pur¬ 
poses we limit ourselves to systems without back mutations. The model in its full generality is 
free of this assumption, although it is unclear how to define pathways when back mutations are 
allowed (see Supplementary Information for a detailed explanation). To introduce the framework 
we begin with a simple model in which the wild type can have two independent mutations leading 
to the fittest type. Then we increase the number of mutational events it takes to get to the corre¬ 
sponding type leading to a generalization of the methodology. We briefly mention an application 
of this approach by linking it to a cancer initiation model [28] showing how mutational epistasis 
changes the path probabilities. Finally we provide an outline on how to extend the model to a 
general system where different mutations need to be acquired to reach the final mutant. 


METHODS AND RESULTS 

Probability Generating Function 

For our methodology, we are making use of extinction probabilities, more specifically the 
probability for different types to be present or not to be present. In a branching process this 
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probability can be recursively obtained using probability generating functions (PGF). Since the 
relation between PGFs and the probability for a type to be present is the main tool we are using, 
we devote this subsection to giving a short overview about this correlation, although it is rather 
technical and well known (e.g. [27, 30]). 

The probability generating function (PGF) in discrete time for a one-type process is in general 
defined as 

OO 

fis) = '^Pks\ ( 1 ) 

k=0 

where k denotes the number of offspring and pk represents the probability of having k offspring 
(the focal individual dies in this context) [27], For many biological processes, for example 
cell multiplication, it makes sense to only consider offspring numbers of 0 (death), 1 (nothing 
happens), and 2 (cell division). But in other biological systems it makes sense to consider 
many offspring at once, for example reproduction via numerous seeds in plants. Our analysis 
is not restricted to any particular offspring distribution. However, for the sake of simplicity, we 
restrict our example to the so called binary splitting, i.e. either two or no offspring. The use of the 
argument s is not obvious at this point. If we set s equal to 0, the probability generating function 
reduces to /(O) = po, which is the extinction probability for a population of one individual in one 
time step. Since all individuals behave independently, /(O)^ = Pq is the extinction probability 
for a population of size N in one time step. Now looking at the extinction probability within two 

time steps, we note that with probability p 2 we would have two individuals in the next time step 

originating from one individual. Hence, the extinction probability for a single individual within 
two time steps is, 

Po+P2pI = /(/(O)) = (2) 

and that of population with N individuals is, 

(po + P2Plf = (/(/(O)))" = (r'^HO))" . (3) 

Continuing for further time steps, we see that /°^*^(0) is the extinction probability for the system 
within t time steps. 

As of now we assumed that individuals reproduce clonally i.e. giving rise to the same type. 
Now we continue investigating the extinction probability for a two-type process. Let us think of 
the two types A and B, where an A individual can produce any number of A or B individuals. 
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and respectively for B. Then the general PGFs if the process starts with one type A or one type 
B individual are defined as 

oo oo 

fA{sA,SB) PkA,kB^'As'B, (4) 

kB=0 
oo oo 

fB{sA.SB)=J2J2pLl.AA4‘’: (5) 

^^=0 kB =0 

where Pk^^kB iPkA,kB) ‘^^'^otes the probability of one A (B) individual producing kA A and ks 
B individuals in the next time step. Let us try to recover the extinction probability as for the 
one-type process. If we set both and sb equal to zero and assume that we start with one A 
individual, we obtain a similar result as above for the total extinction probability 

/.4(0.0)=P^„. (6) 

Oftentimes, one is rather interested in the extinction, or non-presence, of just one particular 
type. Let us for example assume we are only interested in the presence of B individuals. The 
probability of having no B individuals in time step 1 is the sum over all probabilities, where no 
B offspring is being produced J2TA=oPk[^o ^ /a(b)( 1,0), starting with one A (B) individual. 
Now looking at the probability of having no B individuals in time step 2, we need to account for 
the probability of having kA A and ks B individuals being produced in the first time step. This 
leads to 

OO OO 

Y. Y = /.4(/.i(1,0),/b(1,0)) =: /°<">(1,0). (7) 

kA=^ kB=^ 

Continuing this procedure and analogous to the one-type process, the probability of having no 
B individual in time t is /^^*^(1,0). 

In a similar fashion this procedure can be extended to a multi-type process with an arbitrary 
number of types. For further information and detailed insights into extinction of branching 
processes we refer to [27] and [30]. 

Two dimensional fitness landscape 

We begin with a minimal fitness landscape. Envision a wildtype ab which can mutate at the 
two loci to A and B, respectively. With both mutations, the system is in the final state of AB. 
In such a system there are two different paths as illustrated in Figure 1. Traditionally, epistatic 
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FIG. 1. Mutational pathways for a system with two loci. There are two different pathways 
to reach the final mutant. Fitness is represented by the size of the circles denoting the types. 
Thus the wildtype ab and Ab have a similar fitness whereas AB has a significantly greater fitness 
compared to the wildtype while aB is much less fit than the wildtype. When all mutation rates 
are the same, the pathway via aB would be not adaptive, since this type has a low fitness. If 
the mutation rate is large enough, especially if (indicated by the thick arrow), this 

pathway becomes accessible. 

models are discussed in terms of different fitness values, whereas the mutation rates stay the 
same [13, 14]. Exemplarily the fitness landscape for a system with sign epistasis is shown in 
Figure 1. In such a system where the mutation rates stay the same, i.e. pa = Ma = Mb' 

it is clear that the path via Ab is the most probable one. However, if the mutation rates change, 
e.g. Pa, also the path via aB can become accessible. Changing mutation rates amounts 

to including epistasis in the mutational landscape in addition to epistasis in the fitness landscape 

[29]. 

For the four types of the above model, we need to consider four different PGFs, one for each 
type 

fabi^'^ab) '^Abi ^aB) ■^Ab') “F ^a6((l fJ'A fJ'B^^ab “F fJ^A^Ab “F fJ^B^aB^) ) 

/'-Afe('5a6) ^aBi ^Ab') “F nA^Ab ~F fJ'S^AB^ ) 

fasi^^ab^ ^Ab^ ^aB) ^Ab') ”F “F , 

fAB{Sab, -^Ab) ^aB, Sab) = dAB + bABS^si (S) 

where bj and di are the birth and death probabilities of type i. The exponent of 2 arises from 
a branching process with binary splitting. The arguments Sab, ■ ■ ■ ,sab correspond to extinction 
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probabilities of the respective type as discussed above. The functions fi correspond to the 
extinction probability of the whole process given that the process starts with a single individual 
of type i. The PGF fi at time t is recursively calculated as 

fl'\sab, SAb, SaB, Sab) = MfiV\fAb^\faB^\fAB^^)- (9) 


Time Distribution 

Using the generating functions we now approach the extinction time distribution of the binary 
branching process. Particularly starting with 1 wild type individual, the probability of having 
no 74i?-individual at time t is {1,1,1,0) =: f{t). Thus the probability of having at least 
1 ^5-individual at time f is 1 — f{t). The probability, that at least 1 A5-individual appears 
exactly at time t is the probability, that there is an AS-individual at t minus the probability that 
there was already one at time t — 1: 

r{t) = (1 - fit)) - (1 - fit - 1)) = fit - 1) - fit). (10) 

Starting with N wild type individuals the probability that there are no ^5-individual at time t 
is then /(t)^. This leads to the time distribution as, 

r(f) = r(f-i)-r(f). (11) 

However, the arising AB should start a lineage that does not die out. Hence we are interested 
in the probability of having a successful A5-individual. To calculate this we use the known 
extinction probability of an 745-individual in place of sab- The probability of an 745-individual 
going extinct is its death probability divided by its birth probability bab ■= dAs/bAB [!]■ The 
modified PGFs for this purpose then read as 

fabi^ab) ^Ab) ^gb) d^b T ^a6((l I^A t^B)^ab ”F t^A^Ab ”F IJ^B^gb) ; 
fAbi^abi ^Abi ^gb) dAb T ^A6((l t^B^^Ab T I^B^Ab) i 

faBi^ab, ^Ab, ^as) = + ^aB((l “ + fJ^A^AB)"^- (12) 

Note, that the PGF for the final mutant type is not necessary anymore. We can now calculate 
the time distribution until the first successful mutant appears the same way as described above. 
Figure 2 shows the perfect agreement between the recursive solution and 5000 simulations. The 
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parameters, specified in the Figure 2's caption, are entirely arbitrarily chosen to reflect an epistatic 
fitness landscape as sketched in Figure 1. The reason we chose a very slightly advantageous 
fitness for the type 746-individuals is solely to stress the fact, that this method holds for any 
fitness values, not only if some are restricted, for example to being neutral. 



FIG. 2. Time distribution of reaching the final mutant for a four type fitness landscape 

as in Fig. 1. Solid line represents the recursive solution and the bars represent 5000 simulations. 
The parameters are: 

Death probabilities: dab — 0.5^ dAb — 0.49995, — 0.25. Birth probabilities are 

1 minus the corresponding death probability. Mutation probabilities are [ib = I^b ~ 
fiA = ‘2^ ’ 10“^, = 0.005. Population size in the beginning: N = 30000. 


For a three-type continuous time branching process, as in A B C, the time dis¬ 
tribution was computed in [32], This was done using the analytical solution of the probability 
generating function for the two-type process A B [33] and the fact, that in continuous time 
mutations follow a Poisson distribution. Adding a second intermediate type, e.g. B 2 , would also 
give such a process but immediately results in unwieldy analytical calculations. 


Path Probabilities 

In the current example there are two possible paths by which the wildtype can reach the final 
mutant AB, either ab ^ Ab ^ AB or ab aB AB. Experimental evidence shows that 
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not all paths are equally probable [15, 34], Beginning with ab then what is the probability of the 
first AB mutant arising via either path and how long does it take for the different pathways? 

The probability, that the first mutant arises exactly at time t via pathway Ab is (derived in 
the SI), 

P4»(i) = /''(i-l)-(/'"'(t))''. (13) 

where is defined in the Supporting Information (SI) and is being computed in a similar 

fashion as f{t). The total probability for this path is then the summation of pAbit) 

OO 

QAb = ^pAb{t)- (14) 

t=l 

Computationally the sum would go up to a tmax, where ~ 1) ~ f^^’^Ktmax) < ^ 

(where usually machine epsilon is chosen as e). The total extinction probability of a multi¬ 
type branching process is determined by the smallest fixed point s* = {s'^b, of 

the probability generating functions f(s*) = s*, where <b is the extinction probability, if the 
process starts with one a6-individual [27], Nevertheless those total extinction probabilities are 
not suitable for the question, via which path the first successful ^5-mutant arises. The problem 
lies in the time; the pathway via Ab for example could have a very low extinction probability 
whereas the pathway via aB might have an extinction probability of 1/2. Intuitively one would 
expect the path via Ab to be more frequent. However, if the path via aB is much faster (e.g. due 
to S> Pb) would actually find that each path happens with probability that approaches 
1/2. Therefore, it is important to do the recursive analysis to include the probability, that a 
successful mutant did not arise through any other path beforehand. 

Figure 3 shows the probability densities for the different pathways of the minimal model. 
Interestingly, the pathway via aB is predominantly prominent in the beginning but overall less 
likely. Hence if experiments are stopped after a short time interval then they might provide 
conclusions which can be upended by looking at the experiments at a later time point. 

Multiple mutations in two dimensions 

In the earlier model the wildtype had two possible mutations a ^ A and b ^ B. It is possible, 
that a to A and bto B are a multi-step process. Hence we can assume that it takes m mutations 
to go from a to A and n to go from b to B. Hence for m = n = 1 we recover the simple model 
as discussed above. The calculation of the time distribution can be directly transferred from the 
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FIG. 3. Probability distribution for the different pathways. Orange represents the pathway 
via aB and blue the pathway via Ab. The bars are the results of simulations, the solid lines depict 
the computed results. In the pie charts the distribution of the pathways are illustrated up to 500 
time steps (shaded area, left pie chart) and up to 5000 time steps (right pie chart). Stopping after 


a few lineages have reached the final mutant might lead to a false distribution: The other pathway 


might just need longer, but have a smaller extinction probability. The parameters are: Death 
probabilities: dab — 0.5, = IjZ^dAh — 0.49995,^^5 = 0.25. Birth probabilities are 1 minus the 

corresponding death probability. Mutation probabilities are jiB — ~ MA = 2 • 10“^, 

I^A ~ 0.005. Initial Population size is = 30000. 

simple model by including all necessary probability generating functions for all available types. 
Increasing the length of the dimensions has a direct impact on the number of paths leading from 
the wildtype to the final mutant. In particular there are N = possible paths. Assuming 

in general m mutations in the A dimension and n in the B dimension we enumerate the paths 
as follows. Path 1 is the path where at first all A mutations and subsequently all B mutations 
happen. Path 2 is the path where all but one A mutations happen first, then one B, then the 
last A, and finally all other B mutations. Figure 4 shows the different paths for a system with 
four mutations for type A and one mutation for type B. Thus calculating the path probability 
for any particular path p now takes the form, 



( 15 ) 


10 








where f{t) is the probability generating function as in Eq. A.2 and is defined analogously to 
Eq. A.9 in the SI 


Or] 


m+n 


= fi 


PO 


J vn ^ j V 


dn 


fo{t-2) f 

1 J m 1 ■■■ 1 J c 


>(t-2) 


(16) 


Here, the probability generating functions with a p index belong to types along the regarded 
path (which in total are m + n + 1 without back mutations, beginning at 0, with which we 
always label the subindex for the wild type). Accordingly, probability generating functions with 
a q index are associated with types, that do not belong to the respective path (which are in 
total m X n). The probability generating function for the final mutant type is again replaced by 
the extinction probability of this type. We use our framework with this extension on the cancer 
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FIG. 4. Exemplary numbering of the different mutational pathways in a system with 
m — A mutations for type A and n — 1 mutation for B. 


initiation model proposed in [28]. Therein a model with several mutational steps to reach state 
A and one mutational step for state B is analyzed (cf. Fig. 4). The direct change in fitness 
for the A mutations is (nearly) zero, and the B mutation alone is even deleterious. However, if 
an individual obtains all A mutations and the B mutation, the fitness is enhanced which in the 
model leads to rapid proliferation. Here, we provide an example on how the path probabilities 
change, when epistasis is not just in the fitness landscape but in the mutational landscape as 


11 




well. Figure 5 compares the path probability distributions with and without epistasis in the 
mutational landscape. The fitness values, the birth and death probabilities respectively, as well 
as the "nonepistatic" mutation probabilities, are the same as in [28]. 


Multi dimensional fitness landscapes 

The cancer landscape discussed above is a two dimensional system. In principle it is possible 
to extend this approach to higher dimensions. For fitness landscapes of higher orders [15, 35] it is 
still possible to write down the system of probability generating functions and apply the approach 
explained here. The concept remains the same. For each type the probability generating functions 
are needed except for the final mutant type, here only the extinction probability is necessary (SI). 
Finally the probability generating function for the wild type needs to be recursively calculated 
for the time distribution. For the path probabilities the probability generating functions related 
to types not along the considered path again are one time step behind, similar as in Eq. 16. 
However for these experimental fitness landscapes while we can get accurate data elucidating the 
fitness landscape, the mutational landscape is usually hard to determine. 


DISCUSSION 

We have presented a theoretical framework to study mutational pathways in epistatic systems. 
The crucial part is that in our analysis epistasis affects not only fitness (i.e. proliferation and 
death rates) but also mutation rates. Hereby we could show, that pathways become accessible, 
that without mutational epistatic effects are mostly unlikely to happen (cf. e.g. Figure 5). Our 
analysis is based on multi-type branching processes and hence it does not rely on the assumption 
of a constant population size. 

While we have focused on a fairly simple system with a fitness landscape with a single peak, 
the approach can be extended to a rugged fitness landscape. Moreover, if back mutations are 
involved, one can still calculate the time distribution, although pathways are not clearly defined 
in a system with back mutations anymore (see SI). Furthermore in the current scenario in each 
time step the individuals could replicate or die. In addition we could have a resting probability 
where the individuals remain in the same state with a certain probability. Such complicated 
scenarios can be incorporated in our framework as well (SI). The computations can be precisely 
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FIG. 5. Comparison between the path probability distributions of a minimal Burkitt 
Lymphoma model. Top: Time distributions for the model (a) without epistatic effects on 
mutation probabilities and (b) with mutational epistasis. The probability to obtain an A mutation 
is 100 times higher, if the B mutation is present in that individual. Bottom: In (c) the path 
probabilities for the model without epistatic effects on mutations are illustrated, whereas in (d) 
the mutation probability is again increased by 100 for an A mutation if the B mutation is present. 
Pathway 1 corresponds the the mutational pathway, where first all necessary extra mutations have 
to be acquired, and the B mutates last. Pathway 2 denotes the pathway, where 3 of 4 extra 
mutations have been obtained, then the B mutation happens, and at last the final extra mutation 
is acquired. Respectively for the other pathways (cf. Figure 4). The parameters are the same 
as in [28]: The birth probability for an individual with j passenger mutations and without the B 
mutation is 6o,j = 0.5(1 -h 10“^)-^, and with the B mutation • 1.015T The mutation 

probability for the B mutation is fio = ^ ' 10“^, for an A mutation without the B mutation being 
present fip = 2 • 10“^, and with the B mutation being present (only necessary for (b) and (d)) 
fjL^ = 2 • 10“^. The population size in the beginning is N = 500000. 
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represented in analytic terms and need to be solved recursively. 

We apply our framework to a cancer model including mutational epistasis [28] and show how 
the path probabilities are altered by it. Mutational epistasis can thus lead to heterogeneity in 
the density of different mutant types between different age groups as reaching the final mutant 
early is only possible by one mutational pathway which is not possible at later time points. 

As shown here the mutational landscape can undermine the current predictions based solely on 
fitness landscapes. Just like in longterm evolution, experimental as well as theoretical approaches 
ought to be balanced between studying effects of selection and the strengths of mutations. The 
theoretical analysis based on the approach explained here helps in understanding the importance 
of mutational epistasis, even though the computations have to be solved recursively. In particu¬ 
lar, it makes analyzing the fitness and mutational landscapes more interactive, since long-lasting 
simulations are not necessary any more. 
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SUPPORTING INFORMATION: REPEATABILITY OF EVOLUTION ON EPISTATIC 
LANDSCAPES 

General Probability Generating Functions 

In the main text we considered only the case where each individual has to die or divide in 
every time step. Here we relax this assumption and consider a more realistic scenario where only 
some individuals proliferate or die, whereas others do not take any action at all (Fig. 6). Then, 
the probability generating functions for the four types: wild type, individuals with mutation A, 
individuals with mutation B, and individuals with both mutations are defined as 

fabi^abi ^Ab^ ^aBi ^Ab) “h (1- bab + — fJ'A — fJ'B)Sab + t^A^Ab + fJ^B^aB)^, 

fAbi^^abi ^Ab^ ^aBi ’^AB^ ~F (1 -\- 6^^((1 ) 

faBi^ab, SAb) -^aB) -^Ab) = daB + (1 ~ Kb “ daB)SaB + KbH^ ~ l^A)SaB + I^A^AbY) 
fAB^Sabi SAb) ^aB, Sab) = dAB + (1 “ ” dAB)sAB + ^abS^B' ( 1 "^) 

The functions are similar to the scenario of binary splitting (cf. Eq. 1 in the main text). There 

1 — h — d 



FIG. 6. Process described by the general pgf. An individual can either die, proliferate, or 
neither and just live. If it proliferates the offspring can mutate. In case of including back mutations 
additional mutation terms appear leading as in Eq. (18). 

is only one term added: (1 — 6* — di)si,i G {ah, Ah,aB, AB} which denotes the case of the 
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individual neither dividing nor dying. To make the model even more realistic one could also 
include back mutations, 

fab{,^abi ^Abi ^Ab') T (1 ^afe T fJ'A l^B^^ab T fJ'A^Ab T i^B^aB) 

fAb{Sab, SAb, SaB, Sab) = dAb + (1 “ ^Ab “ dAb)sAb + &A6((1 “ l^ab “ I^B)^Ab + l^ab^ab + I^B^AbT 

fasi^^ab) ^Abi ^aBi ^AB^ d^B T (1 ^aB d^B^^aB T ^aB((l l-^ab l^A)^aB T l-^ab^cib ~h l^A^AB^ 
fABi^ab, SAb, SaB, Sab) = dAB + (1 ~ ^ab “ dAB)sAB 

+ bAB ((1 - IJ^A^ - I^B^)saB + l^i^SAb + l^i^SaB)^ 

(18) 

If the fitness landscape is rugged, i.e. having multiple local optima, they would be inaccessible 
from certain “downstream” directions if back mutations are not allowed. Hence allowing back 
mutations, allows to have a rugged fitness landscape with local optima accessible from multiple 
directions. The probability generating functions seem more complex, but the principle of the 
computation as discussed in the main text does not change at all. 


Time distribution 


Here, we give a more detailed description on how to calculate the time distribution for the 
minimal model with four types, and two paths, but with back mutations. 


1. Calculate the extinction probability of the final mutant type AB as in [1] 


dAB + bAB {l^A^ + 


(19) 


Note, that without back mutations the extinction probability reduces to CAB = ^ as in 
the main text. 
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2. Until some t^ax calculate recursively 


,AB ro(t-l) 
laB 


Bab — ^ab + (1 “ ^ab — dAB)fAB 

+ bAB ((1 - 

Bb^ = daB + (1 - baB “ daB)faB~^^ + baB ((1 “ “ ^^aB) faB~^^ + l^aBfab~^^ + /^ a / as ”^^ 

Bm = dAb + (1 - bAb - dAb)BAb~^^ + bAb ((1 - l^aB)BM~^^ + ^^aBBb~^^ + I^bI 

( 20 ) 

fit) := /°f = dab + (1 - bab - dab)f°b~^^ 

+ bab f(l “ /^A — lJ'B)fab + I^a/A b 


A 

AB 


eo{t-l) 

laB 


where = Bb^^ = 1 ^nd = cab- Note, that without back mutations 

these functions would not be coupled anymore and one can first calculate fB and Bb fo*” 
all t, since those functions would not depend on fab- Moreover, would be equal to 
Cab V t. Hence, one would not need to recursively calculate /ab- However, the complexity 
does not change. 


3. The probability to get the final, successful AB mutant, i.e. an individual that produces a 
lineage that does not die out again, exactly at time t is 


r(t) = r(t-i)-r(t). 


( 21 ) 


where N is the number of individuals in the beginning. Calculating this for all t G 
{0,... ,tmax} we obtain the time distribution. 


Single-Path time distribution 

Here, we explain the computation of the probability distribution of the pathway via type Ab 
exemplarily. Allowing back mutations it is unclear how to specify different mutational pathways. 
For instance for the pathway ab ^ aB -A ab ^ Ab ^ AB it is obscure to say via which type 
the final mutant has been reached. Obviously the final mutant has been reached via type Ab, 
but it might be necessary for the population to first reach type aB. Hence, aB might play a 
vital role for reaching AB, too. For this reason we neglect back mutations in the computation 
of the path probabilities, thus guaranteeing clear distinguishable pathways. 
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Let Ab(t) (aB(t)) denote the random variable, that there is an AB mutant until time t via 
pathway Ab (aB). Thus, -iAb(t) corresponds to the random variable, that there is no AB mutant 
until time t vial pathway Ab. Then the probability, that the first mutant arises exactly at time t 
via pathway Ab (i.e. not via pathway aB beforehand) is 

PAb{t) =P{Ab{t) n -<Ab{t — 1) n -<aB{t — 1)) 

=P{^Ab{t - 1) n ^aB{t - 1)) - P{^Ab{t) n ^aB{t - 1)). (22) 

The first term is calculated by the pgf as in Eq. (17). For the second term however, the time 
points for the different pathways are different. Let us derive a recursive function for this second 

term at this point. To do so, let us first consider the extinction probability for the subprocess 

of Ab AB, where the process starts with one Ab individual. As discussed previously, this 
extinction probability within t — 1 time steps can be recursively calculated by its probability 
generating function 

Im + (1 “ ^Ab — dAb)fAb + ^Ab (^(1 “ pB)fAb + Pb^AB^ , (23) 

with — 1. Similarly, the extinction probability for the subprocess aB —> AB within t — 2 
time steps can be calculated recursively using the probability generating function for aB 

faB = daB + (1 “ Kb - daB)faB + ^aB (^(1 “ pT)faB + Pa^AB^ , (24) 

with = 1. When we now consider the extinction probability of the whole process starting 
with an individual of type ab, we see that it can either go extinct right away, or if it divides we 
can refer to the individual extinction probabilities for the different types (in case of mutation), 
i.e. their probability generating functions 

fab^ •= dab + (1 - Kb - dab) fab + Kb (^(1 ” AA “ PB)fab + PaIA b + PsfaB 

_ f (Mt-i) 

— JabKJab ^ J Ab ■> JaB A 

with = 1, = 1, and = 1. Note, that in contrast to the normal probability 

generating function, here the probability generating function for type aB has one time step less, 
which agrees with the second term in 22. To not confuse this modified probability generating 
function with the common one, we use the bar-notation. Again, no probability generating function 
for the AB-type is necessary, since the actual extinction probability for this type is used. 
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We define this recursive function as 

f°h\sab: SAb: SaB, Sab) ■= (26) 

The index Ab denotes, that this is the modified probability generating function for the pathway 
via Ab. 

With this we now describe the algorithm for the path probability. 

1. Calculate the extinction probability of the final mutant type AB as above. 

2. Until some tmax calculate recursively f{t) as explained above in Eq. 20. 

3. Until some tmax calculate recursively 

fa^B = daB + (1 “ Kb - das) faB + ^aB (^(1 ” , 

fAb = dAb + (1 - bAb - dAb)fM~^^ + bAb ((1 - , (27) 

■= IT = dab + (1 - bab - dab)f°T^^ 

+ bab (^(1 - t^A- IAB)fab + I^AfAb + d-sfaB 

where = f~^ = = 1. Note, that the only difference is that the probability 

generating function of types not along the pathway considered is one time step behind 
(marked in red). This is also the reason, why there are two initial conditions needed for 
type aB. 

4. The probability to get the final, successful AB mutant exactly at time t via path Ah and 
not getting a successful AB mutant beforehand is then computed by 

PM, = f(t - 1) - (/''“>(«))" ■ (28) 

Analogously one can calculate the path probability for reaching the final mutant via aB. 
Note, that while this computation gives the correct path probabilities, the sum over all paths 
can be slightly greater than the overall time distribution. This is due to the fact, that in time 
discrete systems the final mutant can be reached by different pathways at the same time. In the 
description here, such cases count for all pathways that succeed at the time. 


[1] K. B. Athreya and P. E. Ney. Branching Processes. Springer, Berlin, 1972. 
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